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Abstract 

Basic ideas of muon tracker technique for the solution of primary cosmic ray (PCR) com¬ 
position problem in the energy range 10^^ — 10^® eV are presented. The approach uses MC 
simulation data made with CORSIKA6.990 for ’’Pamir-XXI” site conditions. Similar technol¬ 
ogy can certainly be developed for other observation levels and interaction models. One can 
probably extend it to much higher primary energies. 


Introduction 

It is useless to repeat all the words that show the importance of the study of primary mass 
composition at super high energies covered by the extensive air shower (EAS) method. Many 
different groups using different detector set-ups have dealt with the problem but the answer is 
still vague. The reason for such a poor result, from our viewpoint, lies in the inadequacy of 
experimental methods used. To be more precise, EAS characteristics that are experimentally 
measured {Xmax, X^j^/Ncharged) do not contain enough information on the mass of the primary 
particle. A natural way out of such a situation is to look for more informative measures of the 
primary mass. 

In 1-100 PeV energy range the EAS Cherenkov light angular distribution characteristics 
prove to be appropriate for the solution of mass composition problem [1,2]. Generally speaking, 
this technique can be used at higher energies as well but the detector array should be made 
much larger while still being rather dense and, thus, must be very expensive. 

Here we present another approach for Eq > 100 PeV which seems to be more cost-effective 
and, probably, will work even for the extremely high energy showers. 


Main idea of the new method 

EAS is a particle cascade originated by a primary cosmic ray, probably proton or other 
nucleus of super high energy. The primary particle undergoes a few strong interactions with air 
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nuclei (N, O, Ar) giving life to quite a number of secondary particles. The most numerous EAS 
component is the electromagnetic one ( 7 ,e''',e“) which is characterized by a radiation length 
(37 g/cm^) that is about 3 times smaller than hadron interaction length. On the contrary, the 
range of the next to most numerous component, which is the muon one, can be much larger 
and for ~ 10 GeV is comparable with the shower length (10-30 km). One can conclude that 
the EAS muons with > 1 GeV can yield important information on the longitudinal shower 
development. In this aspect muons resemble the EAS Gherenkov light which proved to be the 
most informative component of all. Grucial question is how to collect the data on the PGR 
mass composition borne by muons, i.e. what characteristic of the muon component is sensitive 
to the primary mass and, thus, to be measured in EAS experiment. 

Usually, the muon content of the shower is used as a measure of the primary particle mass. 
Unfortunately, any procedure for estimation involves the measurement of muon densities all 
over the area swept by EAS. The measurement is made by local detectors of limited area ~1 
m^ spread widely which results in rather large uncertainties of statistical (due to small detector 
area) and physical (due to cascade fluctuations) origin. Finally, it turns out that A^ uncertainty 
kills its sensitivity to the primary mass. As the ultimate results of KASKADE-Grande [3] show, 
A^ parameter is only capable of dividing all the events into two groups in primary mass which 
is dehnitely not enough to solve the composition problem. 

Gertainly, one should think of Ending more sensitive parameters, preferably directly mea¬ 
surable ones. One of the possibilities is to use the similarity between the muons and Gherenkov 
light. 

An optical detector observing an EAS from large (^300 m) core distance reveals a Gherenkov 
light pulse presenting roughly an electron cascade curve. The latter can also be seen by an 
imaging telescope which field of view is divided into sensitive areas of small size < 1° in diameter 
(pixels). 

Because relativistic muons can reach the observation level even from as high as the first inter¬ 
action of the primary particle, one can think of some muon analogue of the imaging Gherenkov 
telescope, namely, of a muon tracker as a detector of muon angular distribution which should 
be sensitive to the primary mass. A priory it is not clear 

a) whether such a sensor can distinguish even between primary proton and iron nucleus events, 
cascade fluctuations taken into account; 

b) what the design of such a tracker should be; 

c) how must the experimental data on muon angular distribution be processed in order the 
primary mass information to be extracted. 

Full MG simulation presumably can answer the questions. 


Artificial shower generation 

The present study is made in the framework of the ’’Pamir-XXI” project which determines 
the observation level (4250 m above sea level) and the primary energy range of interest (10^® — 
10^® eV for secondary charged particle methods). ’’Pamir-XXI” will include an optical part [2] 
incorporating a net of fast detectors and imaging telescopes that, working together, are capable 
of dealing with PGR mass composition problem in 100 TeV - 100 PeV energy range. At higher 
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energies the event rate with the basic (minimal) optical set-up will be low and one should think 
of a charged particle detector set-up to take over the task. 

In view of all written above, one should first check 10^^ — 10^® eV energy range for a muon 
angular distribution sensitivity to the primary mass. 

The following EAS samples were generated: 

- 100 PeV proton (p), nitrogen (N) and iron nucleus (Fe) initiated vertical showers, sample 
volumes 50; 

- 1 EeV p, N and Fe initiated vertical showers, sample volumes 30. 

For the generation of all samples CORSIKA6.990/QGSJET01 [4] was used without thinning 
with energy threshold 1 GeV for all secondaries. 

Also 100 PeV p, N and Fe initiated vertical shower samples of volume 30 were generated 
by CORSIKA6.990/QGSJET-II to analyze the interaction model effect on muon angular dis¬ 
tribution. 



Figure 1: Detection geometry for vertical shower: polar a and azimuthal 0 of a muon track. 
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Muon angular distribution processing: how to 
distinguish between different primary nuclei 

The following considerations defined the choice of the main parameters of the muon tracker 
and data handling procedure (area, angular resolution, typical core distance, threshold energy, 
azimuthal angle limitations): 

- tracker area should be large enough to detect a number of muons in each shower such that 
a histogram in muon polar angle a (with respect to the shower axis, see Fig.l) can be formed 
with bin contents fluctuating not too much; 


lOOPeV p & Fe vert.EAS at 4250m, muon ADF, R=100m, Ethr=1GeV, phi=-+0.2rad, 50+50 events 



Figure 2: Comparison of muon angular distributions for 100 PeV p and Fe initiated shower 
in case of uniform grid in a (all histogram bins have equal widths). Interaction model: COR- 
SIKA/QGSJETOl . 

- typical core distance R for the primary energies of interest should ensure, on one hand. 
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large enough number of detected muons (i.e. R being not too large), on the other hand, a 
sufficient base to scan the muon cascade curve from the side (i.e. R being not too small); 

- tracker angular resolution should be high enough for the polar angle bins to be made rather 
narrow to scan at least a few kilometers above the observation level, still the tracker must not 
be too expensive due to the angular resolution requirements; 

- energy threshold must ensure that the detected muons can come from high above the obser¬ 
vation level being abundant at the same time; 

- azimuthal angle 0 (see Fig.l) limitations appear due to the fact that the muons of interest 
should come mainly from the shower core and not from its periphery; 

- muon energy threshold and azimuthal angle limitations should not reduce the number of 
detected muons substantially. 

Preliminary tests with ii^o=10 PeV showers showed that the tracker must be of area S ~ 
100 m^ and work at core distance i? > 60 m but the number of detected muons was rather low 
under such conditions {Eq, S, R) to form angular histograms. 

With Eo = 10^^, 10^® eV we use S = 100 m^, R = 100 m within this work but our 
approach can dehnitely be applied to other S and R, the above presented considerations taken 
into account. 

Naive muon angular distribution processing includes histogramming all muons of a sample 
in a above 1 GeV within azimuthal angle 0 band [-0.2,0.2], bins in a being of equal width. An 
example of p and Ee histogram comparison is shown in Fig.2. Marker abscissa correspond to 
the centers of histogram bins, marker ordinates represent bin contents averaged over a sample 
while error bars show r.m.s. deviatons of the contents. One can see no dehnite separation of 
the two distributions. 


r12-r32 vs r15 


r12-r32 vs r14 


r12-r32 vs r13 






Figure 3: Comparison of 100 PeV p and Fe shower samples in different 2-dimensional feature 
spaces in case of uniform grid in a. Features used: and rij — r^j. Interaction model: 

CORSIKA/QGSJETOl . 
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10OPeV p & Fe vert.EAS at 4250m, muon ADF, R=100m, Ethr=1 GeV, opt, 50+50 events 



Figure 4: Comparison of muon angular distributions for 100 PeV p and Fe initiated shower in 
case of optimum grid in a. Interaction model: CORSIKA/QGSJETOl . 


Next step of the naive procedure is the presentation of the two samples in 2-dimensional 
feature space. A few different angular histogram shape measures rij = Wi/wj were tried , here 
Wi is the i-th bin content of a—histogram of a shower. 

The reason for the use of the distribution shape parameters and not absolute muon bin 
contents is simple: relative features are less interaction model dependent and are more robust 
(i.e. less sample dependent). Fig.3 presents some examples of p and Fe samples in 2-dimensional 
feature spaces: one can see substantial sample overlapping in all spaces so that no obvious 
criterion for sample separation can be found. 

Discouraging result of the use of the naive procedure forces one to look for some radical 
improvements of the approach. The first improvement comes from a consideration that equal 
width bins do not reflect the real scale of the muon cascade and one should optimize the bin 
grid. The second one takes into account the steep shape of the a-distribution and rather large 
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p-Fe classification by (w2+w4/w3,w2-w3/w4), 50+50 



p-Fe ciassification by {w2+w4/w3,w2+w3/w4), 50+50 



p-Fe ciassification by (w2-w4/w3,w2+w3/w4), 50+50 



p-Fe ciassification by (w2+w4/w3,w3-w4/w2), 50+50 



p-Fe ciassification by {w2+w4/w3,w2-w4/w3), 50+50 



p-Fe ciassification by {w2+w3/w4,w3-w4/w2), 50+50 



Figure 5: Comparison of 100 PeV p and Fe shower samples in different 2-dimensional feature 
spaces in case of p-Fe-optimum (non-uniform) grid in a. Features used: pf-j^ = ^ 2 +^ _ 
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= P 243 ; etc. Interaction model; CORSIKA/QGSJETOl 


fluctuations of the bin contents: generally speaking, one should reduce the number of histogram 
bins to 4 or 5 leaving the ct-range about the same, which will suppress the fluctuations and 
make the separability of the samples better. 

To carry out the two improvements one is to introduce some quantitative measure of sample 
separability. We use Bhattacharyya distance Db under an assumption that the angular bin 
contents behave as multivariate normal distribution. Each sample is characterized by its own 
mean vector = Wk and covariance matrix k = p,Fe. Then Bhattacharyya distance for 
p-Fe pair takes the form: 


D 


p-Fe 

B 


o (hp k'Fe) S (h; 


k'Fe') T 2^^ 


detTj 


detTjp ■ det'Epe 


S = 


Ej, -|- E 


Fe 


For simplicity one can assume that E^ matrices are diagonal which means the independence of 
different bins of the histogram. The latter assumption is not really true but does not change 
the result substantially. 

After the measure of the class (sample) separability is set one can design a procedure to 
tune the bin width vector h so that Db = Db approaches its maximum for a concrete 

pair of samples. Fig.4 shows the result of such optimization for p — Fe pair of Eq = 10^^ eV. 
One can see dehnite differences between the two histograms as the bin contents (designated by 
markers with error bars) do not overlap for bins 2 to 4, unlike the histograms in Fig.2. Thus, 
we form appropriate histogram shape measures pf-^ = k = 2,3,A,i^j,i^k,j^k 

and present p and Fe samples in different 2-dimensionaI feature spaces in Fig.5. Note that 
the two samples are fully separated for all combinations of pf^j. considered. Fig.6 shows the 
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lOOPeV p, N, Fe vert.EAS at 4250m, muon ADF, R=100m, Ethr=1GeV, opt, 50+50+50 events 



Figure 6: Same as Fig.4 but with N-histogram added. Note that the polar angle grid was 
optimized with respect to p-N separation only. 


same p and Fe histograms as in Fig.4 but with iV-histogram added. The corresponding feature 
space presentations are given in Fig.7. Though the bin widths were optimized for better p — Fe 
separation, A^-sample overlaps p- and Fe-samples only partially for all pf-^ combinations in 
Fig.7. 

One can also optimize the a bin grid for the best separation ofp — iV or N — Fe sample 
pairs. The corresponding feature space plots are presented in Fig.8 and Fig.9 for Eo = 10^^ eV. 
Sample separation is good for all combinations presented. 

We also carry out one more improvement of the bin width optimization algorithm, namely, 
a simultaneous maximization of three pair separation while optimizing the widths. Such im¬ 
provement can be made in many different ways. We do this as follows: 

1. calculate DB,max for each of three pairs p-Fe, p-N, N-Fe D^^rlax) by 

varying the vector h; 










p-N-Fe classification by (w2+w4/w3,w2-w3/w4), 50+50+50 


p-N-Fe classification by (w2-w4/w3,w2+w3/w4), 50+50+50 


p-N-Fe classification by (w2+w4Aiv3,w2-w4/w3], 50+50+50 





p-N-Fe classification by (w2+w4/w3,w2+w3/w4), 50+50+50 p-N-Fe classification by (w2+w4/w3,w3-w4/w2), 50+50+50 p-N-Fe classification by (w2+w3Aiv4.w3-w4/w2), 50+50+50 





Figure 7: Same as Fig.5 but with N-sample added. 


1 EeV p-N classification by (w2+w4Av3,w2-w3/w4). 30+30 


1 EeV p-N classification by (w2-w4Aiv3,w2+w3/w4). 30+30 


1EeV p-N classification by (w2+w4/w3,w2-w4/w3), 30+30 





1 EeV p-N classification by (w2+w4/w3,w2+w3/w4), 30+30 1 EeV p-N classification by (w2+w4/w3,w3-w4/w2), 30+30 1 EeV p-N classification by (w2+w3/w4,w3-w4/w2), 30+30 




Figure 8: Comparison of 1 EeV p and N shower samples in different 2-dimensional feature 
spaces in case of p-N-optimum grid in a. Interaction model: CORSIKA/QGSJETOl . 
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Such algorithm guarantees simultaneous approach oi Dp for each pair to its maximum while 
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1 EeV N-Fe classification by (w2+w4Av3,w2-w3/w4), 30+30 


1 EeV N-Fe classification by (w2-w4Av3,w2+w3/w4), 30+30 


1 EeV N-Fe classification by {w2+w4/w3,w2-w4/w3), 30+30 
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Figure 9: Comparison of 1 EeV N and Fe shower samples in different 2-dimensional feature 
spaces in case of N-Fe-optimum grid in a. Interaction model: CORSIKA/QGSJETOl . 


lOOPeV p-N-Fe separation by (w2+w4/w3,w2-w3/w4), 50+50+50 


lOOPeV p-N-Fe separation by (w2-w4/w3,w2+w3/w4), 50+50+50 


lOOPeV p-N-Fe separation by (w2+w4/w3,w2-w4Av3], 50+50+50 





lOOPeV p-N-Fe separation by (w2+w4/w3,w2+w3Aiv4), 50+50+50 


lOOPeV p-N-Fe separation by (w2+w4/w3,w3-w4/w2), 50+50+50 


tOOPeV p-N-Fe separation by (w2+w3/w4,w3-w4Aiv2), 50+50+50 





Figure 10: Comparison of 100 PeV p, N and Fe shower samples in different 2-dimensional feature 
spaces in case of i?p_Ar_Fe-optimum grid in a. Interaction model: CORSIKA/QGSJETOl . 

keeping the contributions of different pairs in balance. Resulting sample plots for 100 PeV and 
1 EeV showers are presented in Fig.10 and Fig.11, respectively. Sample separation is rather 
good for all 2-dimensional feature spaces considered. 

Finally, we apply the Rp_Ar_i?e-niininiization algorithm to 100 PeV p, N and Fe samples 
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1 EeV p-N-Fe separalion by (w2+w4Aiv3,w2-w3/w4), 30+30+30 


1 EeV p-N-Fe separalion by (w2-w4Aiv3,w2+w3/w4), 30+30+30 


1 EeV p-N-Fe separation by (w2+w4/w3,w2-w4/w3), 30+30+30 



1 EeV p-N-Fe separation by (w2+w4/w3,w2+w3/w4), 30+30+30 


1 EeV p-N-Fe separation by {w2+w4/w3,w3-w4/w2), 30+30+30 


1 EeV p-N-Fe separation by (w2+w3/w4,w3-w4/w2), 30+30+30 



Figure 11: Comparison of 1 EeV p, N and Fe shower samples in different 2-dimensional feature 
spaces in case of i?p_ 7 v-Fe-opfinium grid in a. Interaction model; CORSIKA/QGSJETOl . 


lOOPeV p-N-Fe separation by (w2+w4/w3,w2-w3/wr4), 30+30+30 



lOOPeV p-N-Fe separation by (w2+w4/w3,w2+w3Aiv4], 30+30+30 



lOOPeV p-N-Fe separation by (w2-wr4/w3,w2+w3/wr4), 30+30+30 



lOOPeV p-N-Fe separation by (w2-w3/w1,w3-w4/w1), 30+30+30 



10OPeV p-N-Fe separation by (w2-w3/w1 ,w3-w4/w2), 30+30+30 



lOOPeV p-N-Fe separation by (w2+w3/w4,w3-w4Aiv2), 30+30+30 



Figure 12: Comparison of 100 PeV p, N and Fe shower samples in different 2-dimensional feature 
spaces in case of i?p_ 7 v-Fe-optimum grid in a. Interaction model: CORSIKA/QGSJET-II . 

made with CORSIKA6.990/QGSJET-II in order to check whether the change of the interaction 
model can spoil the muon angular distribution sensitivity to the primary particle mass. Fig. 12 
shows the corresponding sample plots in different 2-dimensional feature spaces. One can see 
that after the change of the interaction model the sensitivity still persists though the set of the 
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most suitable feature spaces may differ for different models. 


Conclusion 

1. A new muon tracker approach to PCR mass composition problem in EAS detection 
is based on high energy (> 1 GeV) muon angular distribution measurements by a large area 
(~100 m^) tracking detector of ~ 0.3° angular resolution placed at core distance R~100 m to 
view the shower sideways. 

2. In order to separate showers initiated by different primary nuclei, a special processing 
of the muon angular distributions is required which maximizes the differences between event 
histograms and, thus, makes it possible to distinguish between different groups of primary 
nuclei. 

3. Concrete realizations of event selection algorithm made for ”Pamir-XXI” conditions 
show that ’’complete” separation of at least three groups of nuclei is possible which makes the 
new method second to none among the charged particle methods in the primary energy range 
considered. 

4. The new approach presumably must work under different observation conditions: obser¬ 
vation levels and/or primary energy ranges. Certainly, the optimal a-histogram bin widths will 
differ in different cases. 

5. The dependence of the method on the hadron/nucleus interaction model should be 
studied next but it is clear now that the sensitivity of muon angular distribution to the primary 
mass will hold because the general properties of the air shower remain the same no matter which 
interaction model is used. 

All in all, large area muon tracker can become an important instrument for the solution of 
the PCR mass composition problem at Eq > 10^^ eV. 
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